# Copyright (c) HySoP 2011-2024
#
# This file is part of HySoP software.
# See "https://particle_methods.gricad-pages.univ-grenoble-alpes.fr/hysop-doc/"
# for further info.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
"""Finite differences stencil inteface used in various numerical methods.
* :class:`~hysop.numerics.stencil.stencil.Stencil`
* :class:`~hysop.numerics.stencil.stencil.CenteredStencil`
* :class:`~hysop.numerics.stencil_generator.StencilGenerator`
* :class:`~hysop.numerics.stencil_generator.CenteredStencilGenerator`
"""
import hashlib
import numpy as np
import scipy as sp
import sympy as sm
import itertools as it
from hysop.tools.htypes import check_instance, first_not_None, to_tuple
from hysop.tools.sympy_utils import recurse_expression_tree
[docs]
class Stencil:
"""
Stencil used for finite abitrary dimension differences.
"""
def __init__(
self,
coeffs,
origin,
order,
dx=sm.Symbol("dx"),
factor=1,
error=None,
delete_zeros=True,
):
"""
Stencil used for finite differences schemes (n-dimensional).
Parameters
----------
coeffs: array_like
Coefficients of the stencil.
origin: array_like
Origin of the stencil from left-most point.
order: int
Order of the stencil.
Other Parameters
----------------
dx: sympy.Symbol or array_like of sympy.Symbol
Spatial spacing symbols.
factor: dtype or sympy.Expr
Common leading factors for all coefficients.
error: sympy.Expr
Exact error if known.
Attributes
----------
coeffs: np.ndarray
Coefficients of the stencil.
origin: np.ndarray
Origin of the stencil from left-most point.
order: int
Spatial order of the stencil.
dim: int
Spatial dimension of the stencil.
shape: np.ndarray
Shape of the stencil.
L: np.ndarray
Distance from origin to left-most point included in stencil support.
R: np.ndarray
Distance from origin to right-most point included in stencil support.
dx: sympy.Symbol or array_like of sympy.Symbol
Spatial spacing symbols.
factor: dtype or sympy.Expr
Leading factor of all coefficients.
error: sympy.Expr
Exact error term if known or else None.
Raises
------
ValueError
Raised if one component of `origin` is negative.
See Also
--------
:class:`CenteredStencil`: Centered version of a stencil.
:class:`StencilGenerator`: Generate Stencil objects.
"""
coeffs = np.atleast_1d(coeffs)
origin = np.atleast_1d(origin)
order = np.atleast_1d(order)
dx = np.atleast_1d(dx)
if (origin < 0).any():
raise ValueError(f"Origin component < 0!\norigin={origin}")
self.dx = dx[0] if dx.size == 1 else dx
self.error = error
self.origin = origin[0] if origin.size == 1 else origin
self.order = order[0] if order.size == 1 else order
self.factor = factor
if delete_zeros:
coeffs = self._delete_zeros(coeffs)
self.coeffs = coeffs
self._update_attributes()
[docs]
def reshape(self, new_shape):
"""Reshape a stencil by adding zeros."""
new_shape = np.asarray(new_shape)
shape = np.asarray(self.shape)
assert new_shape.ndim == shape.ndim
assert (new_shape >= shape).all()
assert ((new_shape - shape) % 2 == 0).all()
zeros = (new_shape - shape) // 2
slc = tuple(slice(z, z + s, 1) for (z, s) in zip(zeros, shape))
new_origin = zeros + self.origin
new_coeffs = np.zeros(shape=new_shape, dtype=self.coeffs.dtype)
new_coeffs[slc] = self.coeffs
return self.__class__(
coeffs=new_coeffs,
origin=new_origin,
order=self.order,
dx=self.dx,
factor=self.factor,
error=self.error,
delete_zeros=False,
)
[docs]
def has_factor(self):
return self.factor != 1
[docs]
def non_zero_coefficients(self):
return np.sum(self.coeffs != 0)
[docs]
def replace_symbols(self, dic):
if isinstance(self.factor, sm.Basic):
self.factor = self.factor.xreplace(dic)
coeffs = self.coeffs.ravel()
for i, coeff in enumerate(coeffs):
if isinstance(coeff, sm.Basic):
coeffs[i] = coeff.xreplace(dic)
[docs]
def apply(self, a, out=None, axis=None, symbols=None, iview=Ellipsis):
"""
Apply this stencil the input array a.
Boundaries are considered periodic.
"""
symbols = first_not_None(symbols, {})
check_instance(a, np.ndarray)
check_instance(symbols, dict)
adim = a.ndim
sdim = self.dim
assert sdim <= adim, "Stencil dimension greater than array dimension."
assert (
set(symbols.keys()) == self.variables()
), f"Missing symbols {self.variables()-set(symbols.keys())}."
out = first_not_None(out, np.empty_like(a[iview]))
axis = first_not_None(to_tuple(axis), tuple(range(adim))[-sdim:])
assert len(axis) == sdim
assert out.ndim == a.ndim
assert out.shape == a[iview].shape
out[...] = 0
for shift, coeff in self.items(include_factor=False):
if isinstance(coeff, sm.Basic):
coeff = coeff.xreplace(symbols)
coeff = float(coeff)
view = a.view()
for i, si in enumerate(shift):
view = np.roll(view, shift=-si, axis=axis[i])
out[...] += coeff * view[iview]
factor = self.factor
if isinstance(self.factor, sm.Basic):
factor = self.factor.xreplace(symbols)
factor = float(factor)
out[...] *= factor
return out
[docs]
def __call__(self, *args, **kwds):
"""Alias for apply."""
return self.apply(*args, **kwds)
def _update_attributes(self):
self.dim = self.coeffs.ndim
self.shape = self.coeffs.shape
self.L = np.asarray([self.origin] if np.isscalar(self.origin) else self.origin)
self.R = self.shape - self.origin - 1
def _delete_zeros(self, coeffs):
dim = coeffs.ndim
shape = list(coeffs.shape)
dtype = coeffs.dtype
if dtype in [np.float16, np.float32, np.float64]:
mask = np.isclose(coeffs, 0, atol=2 * np.finfo(dtype).eps)
else:
mask = coeffs == 0
keep_mask = np.ones_like(coeffs, dtype=bool)
for d in range(dim):
laccess = [slice(None)] * dim
raccess = [slice(None)] * dim
laccess[d] = slice(0, 1)
raccess[d] = slice(-1, None)
laccess = tuple(laccess)
raccess = tuple(raccess)
ldel = mask[laccess].all()
rdel = mask[raccess].all()
if ldel:
keep_mask[tuple(laccess)] = False
if dim == 1:
self.origin -= 1
else:
self.origin[d] -= 1
shape[d] -= 1
if rdel:
shape[d] -= 1
keep_mask[tuple(raccess)] = False
coeffs = coeffs[keep_mask].reshape(shape)
return coeffs
[docs]
def accuracy(self):
"""
Accuracy of the stencil with bigO notation.
Returns
-------
accuracy: sympy.Expr
Accuracy in bigO notation.
"""
expr = 0
done = []
if self.dim != 1:
for dx, order in zip(self.dx, self.order):
if dx not in done:
expr += sm.O(dx ** int(order))
done.append(dx)
else:
expr += sm.O(self.dx ** int(self.order))
return expr
[docs]
def variables(self):
"""
Symbolic variables used by the stencil.
Returns
-------
vars: set of sympy.Expr
"""
_vars = set()
def op(expr):
if isinstance(expr, sm.Symbol):
_vars.add(expr)
recurse_expression_tree(op, self.factor)
for coeff in self.coeffs.flatten():
recurse_expression_tree(op, coeff)
return _vars
[docs]
def coo_stencil(self):
"""
Return a 2d stencil as a sparse coordinate matrix.
Returns
-------
coo_stencil: scipy.sparse.coo_matrix
Sparse Coordinate Matrix of the stencil.
Raises
------
RuntimeError
Raised if stencil is not 2d.
"""
if self.dim != 2:
raise RuntimeError("Stencil is not 2d !")
return sp.sparse.coo_matrix(self.coeffs, shape=self.shape, dtype=self.dtype)
[docs]
def items(self, include_factor=True):
"""
Return an (offset,coefficient) iterator iterating on all **non zero** coefficients.
Offset is taken from origin.
Returns
-------
it : itertools.iterator
Zipped offset and coefficient iterator.
"""
factor = self.factor if include_factor else 1
def mapfun(x):
offset = x - self.origin
value = factor * self.coeffs[x]
return (offset, value)
iterator = np.ndindex(self.shape)
iterator = map(mapfun, iterator)
iterator = filter(lambda x: x[1] != 0, iterator)
return iterator
[docs]
def refactor(self, factor):
"""
Factorize the stencil by a new numeric or symbolic value.
Parameters
----------
factor : :attr:`dtype` or sympy.Expr
New factor to be applied.
Returns
-------
self : :class:`Stencil`
"""
old_factor = self.factor
cfactor = old_factor / factor
coeffs = self.coeffs
for idx, coeff in np.ndenumerate(coeffs):
coeffs[idx] = sm.simplify(coeff * cfactor)
self.factor = factor
self.coeffs = coeffs
return self
[docs]
def hash(self, keep_only=16):
"""
Hash the stencil with its origin, order and coefficients.
The hash algorithm used is sha1.
By default only the 16 first chars are kept from the generated hash.
Parameters
----------
keep_only : int
Number of chars kept from the hashed hexedecimal string.
Returns
-------
hash : string
Hexadecimal hash string of size :attr:`keep_only`.
"""
coeffs = self.coeffs
m = hashlib.sha1()
m.update(self.origin)
m.update(self.order)
for d in coeffs:
m.update(str(d).encode("utf-8"))
m.update(str(self.dx).encode("utf-8"))
return m.hexdigest()[:keep_only]
[docs]
def is_symbolic(self):
"""
Check if either the stencil leading factor or the stencil coefficients
contains symbolic values.
Returns
-------
is_symbolic : bool
"""
return len(self.variables()) != 0
[docs]
def is_centered(self, axe=None):
"""
Check if the stencil is centered (ie. stencil generate a centered
finite differences scheme).
Returns
-------
is_centered : bool
"""
L, R = self.L, self.R
if axe == None:
return (L == R).all()
else:
return (L[axe] == R[axe]).all()
[docs]
def is_cross(self):
"""
Check if the stencil is a cross (zeros everywhere except on a
n-dimensional axis aligned cross).
Examples:
oo+oo o+oo
oo+oo ++++
+++++ o+oo
oo+oo o+oo
Returns
-------
is_cross : bool
"""
mask = np.ones(self.shape, dtype=bool)
if self.dim == 1:
return True
else:
access = self.origin.tolist()
for i in range(self.dim):
acc = [x for x in access]
acc[i] = slice(0, self.shape[i])
mask[acc] = False
return not np.any((self.coeffs != 0) * mask)
def __str__(self):
return """
== {}D Stencil ==
origin: {}
order: {}
accuracy: {}
error: {}
shape: {}
centered: {}
cross: {}
symbolic: {}
variables: {}
factor: {}
coeffs:
{}
=================
""".format(
self.dim,
self.origin,
self.order,
self.accuracy(),
self.error,
self.shape,
self.is_centered(),
self.is_cross(),
self.is_symbolic(),
self.variables(),
self.factor,
self.coeffs,
)
[docs]
class CenteredStencil(Stencil):
"""
Centered stencil used for finite abitrary dimension differences.
"""
def __init__(
self, coeffs, origin, order, dx=sm.Symbol("dx"), factor=1, error=None, **kwds
):
"""
Centered stencil used for finite abitrary dimension differences.
See also
--------
:class:`CenteredStencilGenerator`: Centered stencil generator.
"""
shape = np.asarray(coeffs.shape)
if (shape % 2 == 0).any():
raise ValueError("Shape compnonent even!")
if (origin != (shape - 1) // 2).any():
print(origin)
print((shape - 1) // 2)
raise ValueError("Origin is not centered!")
super().__init__(coeffs, origin, order, dx, factor, error, **kwds)
[docs]
def is_centered(self):
return True
[docs]
@staticmethod
def from_stencil(stencil):
s = stencil
if not s.is_centered():
raise ValueError("Given stencil is not centered!")
return CenteredStencil(s.coeffs, s.origin, s.order, s.dx, s.factor, s.error)